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Summary 


The  overall  objective  of  the  research  carried  out  over  last  year  was  the  development  of 
new  algorithms  for  the  efficient  simulation  of  viscous  compressible  flows  with  moving 
structurally  and  thermally  responding  bodies  in  three  dimensions  using  unstructured 
grids.  The  development  was  based  on  current  3-D  Euler /Navier-Stokes  capabilities, 
and  encompassed: 

-  Grid  Generation; 

-  Flow  Solvers; 

-  Adaptive  embedded  unstructured  grid  methods; 

-  Optimal  shape  design. 

These  topics  are  treated  below  in  some  more  detail. 

1.  GRID  GENERATION 

Many  simulation  techniques  in  computational  mechanics  require  a  space-filling  cloud  of 
arbitrary  objects.  For  the  case  of  ‘gridless’  or  ‘mesh  free’  partial  differential  equation 
(PDE)  solvers,  these  are  simply  points.  For  discrete  element  methods  these  could  be 
spheres,  ellipsoids,  polyhedra,  or  any  other  arbitrary  shape.  The  task  is  therefore  to 
fill  a  prescribed  volume  with  these  objects  so  that  they  are  close  but  do  not  overlap  in 
an  automatic  way. 

Several  techniques  have  been  used  to  place  these  objects  in  space.  The  so-called‘fill  and 
expand’  or  ‘popcorn’  technique  starts  by  first  generating  a  coarse  mesh  for  the  volume  to 
be  filled.  This  subdivision  of  the  volume  into  large,  simple  polyhedra  (elements),  is,  in 
most  cases  performed  with  hexahedra.  The  objects  required  (points,  spheres,  ellipsoids, 
polyhedra,  etc.)  are  then  placed  randomly  in  each  of  these  elements.  These  are  then 
expanded  in  size  until  contact  occurs  or  the  desired  fill-ratio  has  been  achieved.  An 
obvious  drawback  of  this  technique  is  the  requirement  of  a  mesh  generator  to  initiate 
the  process.  A  second  class  of  techniques  are  the  ‘advancing  front’  or  ‘depositional’ 
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methods.  Starting  from  the  surface,  objects  are  added  where  empty  space  still  exists.  In 
contrast  to  the  ‘fill  and  expand’  procedures,  the  objects  are  packed  as  close  as  required 
during  introduction.  Depending  on  how  the  objects  are  introduced,  one  can  mimic 
gravitational  or  magnetic  deposition,  layer  growing,  or  size-basedgrowth.  Furthermore, 
so-called  radius  growing  can  be  achieved  by  first  generating  a  coarse  cloud  of  objects, 
and  then  growing  more  objects  around  each  of  these.  In  this  way,  one  can  simulate 
granules  or  stone. 

We  developed  a  scheme  that  allows  for  the  direct  generation  of  clouds  of  arbitrary 
objects  with  the  same  degree  of  flexibility  as  advanced  unstructured  mesh  generators. 
The  mean  distance  between  objects  (or,  equivalently,  the  material  density)  is  specified 
by  means  of  background  grids,  sources  and  densityattached  to  CAD-entities.  In  order 
not  to  generate  objects  outside  the  computational  domain,  we  assume  an  initial  tri¬ 
angulation  of  the  surface  that  is  compatible  with  the  desired  mean  distance  between 
objects  specified  by  the  user.  Starting  from  this  initial  ‘front’  of  objects,  new  objects 
are  added,  until  no  further  objects  can  be  introduced.  Whereas  the  advancing  front 
technique  for  the  generation  of  volume  grids  removes  one  face  at  a  time  to  generate 
elements,  the  present  scheme  removes  one  object  at  a  time,  attempting  to  introduce  as 
many  objects  as  possible  in  its  immediate  neighbourhood. 

2.  FLOW  SOLVERS 

Any  coupled  fluid-structure  capability  requires  an  advanced,  robust  flow  solver  that 
works  optimally  across  the  range  of  physical  phenomena  and  associated  spatial  and 
temporal  scales.  We  have  made  improvements  on  simulating  three-dimensional  com¬ 
pressible  turbulent  flows  for  complex  geometries. 

Computational  fluid  dynamics  (CFD)  has  become  an  indispensable  tool  for  the  aerody¬ 
namic  design  and  analysis  of  aircrafts  and  engines.  The  difficulty  of  generating  struc¬ 
tured  grids  and  the  desire  to  compute  flows  over  increasingly  complex  configurations 
has  fueled  interest  in  the  development  of  unstructured  grid  methods.  Unstructured 
grids  provide  great  flexibility  in  tackling  complex  geometries  encountered  in  practice 
and  offer  a  natural  framework  for  adapting  to  flow  features.  Nowadays,  unstructured 
grids  composed  of  geometric  simplices,  i.e.,  triangle  in  two  dimensions  and  tetrahedra 
in  three  dimensions  are  widely  and  routinely  used  for  producing  high-quality  solutions 
for  both  inviscid  flows  and  low  Reynolds  number  viscous  flows.  However,  so  far,  the 
solution  of  the  Navier-Stokes  equations  for  high  Reynolds  number  viscous  flows  has 
been  dominated  by  structured  grid  methods,  as  grid  stretching  along  shear  layers  is  a 
trivial  task  on  structured  meshes  but  is  relatively  difficult  to  achieve  on  unstructured 
triangular  or  tetrahedral  grids.  One  popular  approach  has  been  to  use  a  hybrid  grid: 
prismatic  elements  are  used  for  viscous  regions  where  the  grid  is  highly  stretched,  and 
tetrahedral  elements  for  inviscid  regions  where  the  grid  is  nearly  isotropic.  Although 
this  approach  seems  to  be  attractive,  the  need  to  generate  a  prismatic  grid  unavoidably 
renders  the  grid  generation  more  difficult  and  less  robust  for  complex  geometries,  which 
would  otherwise  has  been  the  main  advantage  of  unstructured  grid  methods.  In  fact, 
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one  can  easily  show  that  the  prismatic  grid  generator  will  fail  for  the  corners,  where 
the  newly  created  point  is  not  visible  from  all  the  faces  surrounding  the  original  point. 
Furthermore,  the  existence  of  several  cell  types  in  a  hybrid  grid  undoubtedly  makes 
mesh  generation  and  adaptation,  flow  solution,  and  post-processing  more  complex. 

We  have  developed  a  parallel,  accurate,  fast,  matrix-free  implicit  method  to  solve  com¬ 
pressible  turbulent  flow  problems  using  the  Spalart  and  Allmaras  one  equation  turbu¬ 
lence  model  on  unstructured  meshes.  In  this  method,  the  mean-flow  and  turbulence- 
model  equations  are  decoupled  in  the  time  integration  in  order  to  facilitate  the  incor¬ 
poration  of  different  turbulence  models  and  reduce  memory  requirements.  Both  mean 
flow  and  turbulent  equations  are  integrated  in  time  using  a  linearized  implicit  scheme. 
A  fast,  matrix-free  implicit  method,  GMRES-fLU-SGS,  is  then  applied  to  solve  the 
resultant  system  of  linear  equations.  The  spatial  discretization  is  carried  out  using  a 
hybrid  finite  volume  and  finite  element  method,  where  the  finite  volume  approximation 
based  on  a  containment  dual  control  volume  rather  than  on  the  more  popular  median- 
dual  control  volume  is  used  to  discretize  the  inviscid  fluxes,  and  the  finite  element 
approximation  is  used  to  evaluate  the  viscous  flux  terms. 

The  developed  method  has  been  used  to  compute  a  variety  of  compressible  turbulent 
flow  problems  over  a  wide  range  of  flow  conditions.  The  computational  results  demon¬ 
strated  that  this  unstructured  tetrahedral  grid  method  is  able  to  produce  same  quality 
solutions  as  its  hybrid  conuterpart. 

The  developed  method  has  been  used  to  predict  drags  in  the  transonic  regime  for  both 
DLR-F4  and  DLR-F6  configurations  to  assess  the  accuracy  and  efficiency  of  the  method. 
The  results  obtained  are  in  good  agreement  with  experimental  data,  indicating  that 
the  present  method  provides  an  accurate,  efficient,  and  robust  algorithm  for  computing 
turbulent  flows  for  complex  geometries  on  unstructured  grids. 

3.  ADAPTIVE  EMBEDDED  UNSTRUCTURED  GRID  METHODS 

The  numerical  solution  of  Partial  Differential  Equations  (PDEs)  is  usually  accomplished 
by  performing  a  spatial  and  temporal  discretization  with  subsequent  solution  of  a  large 
algebraic  system  of  equations.  The  spatial  discretization  is  commonly  performed  via 
polyhedra,  also  called  (finite)  volumes  or  elements.  The  final  assembly  of  these  polyhe- 
dra  yields  the  so-called  mesh.  The  transition  from  an  arbitrary  surface  description  to  a 
proper  mesh  still  represents  a  difficult  task.  Considering  the  rapid  advance  of  computer 
power,  together  with  the  perceived  maturity  of  field  solvers,  an  automatic  transition 
from  arbitrary  (and  possibly  imperfect  or  dirty)  surface  description  to  mesh  becomes 
mandatory. 

Two  types  of  grids  are  most  commonly  used:  body-conforming  and  embedded.  For 
body-conforming  grids  the  external  mesh  faces  match  up  with  the  surface  (body  sur¬ 
faces,  external  surfaces,  etc.)  of  the  domain.  This  is  not  the  case  for  the  embedded 
approach  (also  known  as  ficticious  domain,  immersed  boundary  or  Cartesian  method), 
where  the  surface  is  placed  inside  a  large  mesh  (typically  a  regular  parallelepiped),  with 
special  treatment  of  the  elements  close  to  the  surfaces. 
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Prompted  by  the  need  to  solve  shock-structure  interaction  problems  and  the  inability 
of  Computational  Structural  Dynamics  (CSD)  codes  to  ensure  strict  no-penetration 
during  contact,  a  simple  embedded  technique  operating  on  adaptive,  unstructured  grids 
was  developed  in  FY  2003.  The  essential  elements  of  this  technique  may  be  summarized 
as  follows: 

a)  The  key  modification  of  the  original,  body  fitted  edge-based  solver  was  the  removal 
of  all  geometry-parameters  (essentially  the  area  normals)  belonging  to  edges  cut  by 
embedded  surface  faces. 

b)  Several  techniques  to  improve  the  treatment  of  boundary  points  close  to  the 
immersed  surfaces  were  implemented.  Higher-order  boundary  conditions  were  also 
achieved  by  duplicating  crossed  edges  and  their  endpoints. 

c)  Geometric  resolution  and  solution  accuracy  were  enhanced  by  adaptive  mesh  re¬ 
finement  that  was  based  on  the  proximity  to  or  the  curvature  of  the  embedded  CSD 
surfaces. 

d)  In  order  to  save  work,  user-defined  or  automatic  deactivation  for  regions  inside 
immersed  solid  bodies  was  employed. 

The  technique  has  been  extremely  successful  for  complex  fluid-structure  interaction 
cases,  as  well  as  some  external  aerodynamics  cases.  Led  on  by  this  success,  the  tech¬ 
nique  has  seen  extensive  use,  leading  to  new  developments.  These  developments  may 
be  grouped  as  follows: 

-  Improvements  in  speed,  with  particular  emphasis  on  bin  data  structures; 

-  Treatment  of  multimaterial  problems,  in  particular  water/air; 

-  Particle  transport  and  embedded  surfaces; 

-  Direct  link  to  Discrete  Particle  Methods  (DPM); 

-  Body-fitted  grid  generation  via  the  embedded  approach;  and 

-  Links  to  simplified  CSD  models. 

4.  OPTIMAL  SHAPE  DESIGN 

The  availability  of  accurate,  reliable  and  fast  CFD  solvers  immediately  leads  to  the 
next  stage  in  simulations:  optimal  shape  design.  It  is  the  nature  of  any  technical 
development  that  an  engineer  would  not  only  like  to  know  if  a  certain  device  performs 
as  required,  but  how  to  improve  it.  One  of  the  most  promising  ways  of  solving  the 
optimization  problem  is  via  the  so-called  adjoint  formulation.  We  have  developed  a 
scheme  based  on  the  solution  of  a  continuous  adjoint  problem.  The  computation  of 
the  objective  function  sensitivities  is  performed  by  using  only  the  flow  and  adjoint 
solutions  on  the  boundary.  This  property  makes  the  scheme  totally  independent  of 
the  numerical  solvers  (flow  and  adjoint).  In  addition,  the  scheme  allows  the  use  of 
incomplete  gradients  by  only  taking  into  account  the  geometrical  sensitivities,  and/or 
the  boundary  conditions  of  the  adjoint  problem.  An  improved  pseudo-shell  approach 
to  move  the  surface  to  be  optimized,  and  to  impose  the  geometric  restrictions  was 
developed  and  tested.  We  have  incorporated  above  adjoint  solvers  into  FEFLO. 
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The  new  developed  adjoint  solver  was  used  for  an  incompressible  Navier-Stokes  control 
problem,  and  a  compressible  shape  design  problem.  The  control  problem  consists  of 
minimizing  the  amplitude  of  the  periodic  pressure  lift  force  on  the  cylinder,  by  moving 
a  portion  of  the  cylinder  surface  in  each  time  step,  or,  instead  of  that,  by  imposing 
an  equivalent  injection  velocity.  The  control  starts  to  be  applied  at  t  =  25  (500  time 
steps).  The  amplitude  of  the  force  is  reduced  by  87%  for  both  types  of  controls.  In 
addition,  it  can  be  noticed  that  the  effect  of  moving  the  mesh  and  of  applying  the 
injection  velocity  is  exactly  the  same.  The  maximum  injection  velocity  was  around  1% 
of  the  inflow  velocity. 

To  illustrate  how  the  pseudo-shell  approach  works,  the  drag  force  at  constant  lift  over 
a  hypersonic  wing  (at  an  angle  of  attack  of  5°)  was  optimized,  subject  to  a  minimum 
thickness  constraint.  It  was  demonstrated  from  this  numerical  example  that  an  im¬ 
proved  pseudo-shell  approach,  to  move  the  objective  surface,  and  to  impose  general 
geometric  constraints  produced  smooth  shapes  in  each  design  cycle. 


5.  EXAMPLES 

The  developed  unstructured  grid  solver  has  applied  over  the  last  year  to  a  number  of 
problems.  We  include  here  two  recent  examples  in  the  following  section. 

5.1  DLR-F4  configuration:  The  first  example  considers  the  test  cases  for  the  first  AIAA 
Drag  Prediction  Workshop  held  in  Anaheim,  CA  June  2001.  The  workshop  challenge 
was  to  compute  the  lift,  drag,  and  pitching  moment  for  the  DLR-F4  wing-body  con¬ 
figuration  for  different  sets  of  conditions.  For  all  cases  the  free  stream  Mach  number 
is  0.75,  and  the  Reynolds  number  is  3  millions.  The  mandatory  calculations  for  the 
Drag  Prediction  Workshop  were  performed  using  the  standard  unstructured  nodal  grid, 
which  contains  9,650,674  elements,  1,641,451  points,  and  48,339  boundary  points,  as 
shown  in  Figure  la.  All  computations  were  initiated  as  a  uniform  flow  at  freestream 
conditions,  and  advanced  with  a  CFL  number  of  200.  The  relative  L2  norm  of  the  den¬ 
sity  residual  is  taken  as  the  criterion  for  convergence  history.  The  solution  tolerance  for 
GMRES  was  set  to  0.1  with  10  search  directions  and  20  iterations.  We  observed  that 
during  the  first  few  time  steps,  more  iterations  are  required  to  solve  the  system  of  the 
linear  equations:  even  20  iterations  can  not  guarantee  that  the  stopping  criterion  will 
be  satisfied  for  some  cases.  However,  it  takes  less  than  20  iterations  to  solve  the  linear 
equations  at  a  later  time,  and  global  convergence  is  not  affected  by  a  lack  of  linear 
system  convergence  during  the  first  few  time  steps.  All  computations  were  performed 
on  a  SGI  origin  2000  (400MHz)  computer  and  required  4  GB  memory  and  about  7 
hours  on  16  processors  for  1000  time-steps.  Note  that  convergence  in  the  lift  and  drag 
is  obtained  in  as  little  as  500  time  steps  for  all  cases. 

The  case  with  a  single  point  at  Cl— 0.5  is  presented  here.  The  computed  pressure 
contours  displayed  in  Fig.  lb,  clearly  show  good  resolution  of  the  upper  surface  shock. 
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Fig.  la:  Unstructured  surface  mesh  used  for  computing  turbulent  flow  past  DLR-F4 
configuration  (nelem— 9, 650,674,  npoin— 1,641,451,  nboun— 48,339). 


Fig.  lb:  Computed  pressure  contours  on  the  surface  of  DLR-F4  configuration  at 
=  0.75,  a  =  -0.002°,  and  Re-3,000,000. 

Table  1  summarizes  the  results  obtained  for  case  1  by  the  experiments,  drag  predic- 
tion  workshop  participants,  and  present  computation.  Although  there  is  a  considerable 
amount  of  scatter  in  the  data  obtained  by  the  drag  prediction  workshop,  our  compu¬ 
tational  results  are  close  to  the  average  data,  especially  for  drag  and  pitching  moment. 


Table  1.  Comparison  of  computational  results  with  experimental  data 
and  Drag  Prediction  Workshop  data 


Data 

a 

CL 

Cd 

Cm 

ONERA 

0.192 

0.50 

0.02896 

-0.1260 

NLR 

0.153 

0.50 

0.02889 

-0.1301 

DRA 

0.179 

0.50 

0.02793 

-0.1371 

Workshop(Ave) 

-0.237 

0.5002 

0.03037 

-0.1559 

Workshop(Min) 

-1.000 

0.4980 

0.02257 

-0.2276 

Workshop(Max) 

1.223 

0.5060 

0.04998 

0.0481 

Present 

-0.002 

0.5009 

0.03058 

-0.1513 

5.2  DLR-F6  configuration:  The  second  example  considers  the  test  cases  in  the  2nd 
AIAA  Drag  Prediction  Workshop  held  in  Orlando,  FL  June  2003.  The  chosen  con¬ 
figuration,  denoted  as  DLR-F6,  consists  of  wing-body  geometry  without  or  with  a 
conventional  nacelle.  All  computations  were  performed  on  the  NAVO  IBM  Cluster 
1600  (1.3  GHz  Power  4)  computer  and  required  less  than  8  Gbytes  memory  and  less 
than  12  hours  of  CPU  on  8  processors.  Although  the  Drag  Prediction  Workshop  speci¬ 
fied  several  mandatory  and  optional  cases,  only  the  mandatory  cases  were  computed  in 
this  work.  For  all  cases  the  free  stream  Mach  number  is  0.75,  and  the  Reynolds  number 
based  on  reference  chord  of  14.12cm  is  3  millions. 

The  case  with  a  single  point  at  a  prescribed  lift  coefficient  of  0.5.  is  presented  here.  In 
order  to  test  the  robustness  of  the  present  method  on  grids  generated  by  different  grid 
generators,  the  grid  convergence  study  was  performed  using  the  coarse  grid  provided  by 
NASA  LARC  as  our  coarse  grid,  the  coarse  grid  provided  by  Swansea  as  our  medium 
grid,  and  the  coarse  grid  by  DLR  as  our  fine  grid.  Table  2  shows  the  grid  size  and 
computed  results  for  the  wing-body  configuration  for  case  1  with  the  three  different 
grids.  The  comparison  of  lift  and  drag  coefficients  with  the  experimental  data  is  shown 
in  Figs.  2a-b,  respectively.  Although  there  is  a  considerable  amount  of  scatter  in  the 
data  obtained  by  the  drag  prediction  workshop,  our  computational  results  are  close  to 
the  average  data.  The  computed  pressure  contours  on  the  medium  mesh  is  displayed 
in  Fig.  2c. 

Table  2  shows  the  grid  size  and  computed  results  for  the  wing-body-pylon-nacelle  con¬ 
figuration  with  three  different  grids. 


Table  2.  Wing/body/pylon/nacelle  Configuration 


Grid  Size 

a 

CL 

CD 

Cm 

1,827,470 

0.950 

0.497 

0.03719 

-0.1205 

2,419,388 

0.843 

0.498 

0.03658 

-0.1216 

3,682,535 

0.800 

0.501 

0.03631 

-0.1242 
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Lift  coefficient  comparison  with  experimental  data 
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Fig.  2a:  Comparison  of  computed  lift  coefficient  and  experimental  data  for  DLR-F6 
wing-body  configuration  at  Moo  =  0.75,  and  Re=3,000,000. 


Drag  polar  comparison  with  experimental  data 
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Fig.  2b:  Comparison  of  computed  and  experimental  drag  polar  for  DLR-F6  wing-body 
configuration  at  =  0.75,  and  Re=3,000,000. 
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Fig.  2c:  Computed  pressure  contours  on  the  surface  of  DLR-F6  wing-body  configuration 
at  Mqo  =  0.75,  Cl=0.5,  and  Re=3,000,000. 
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5.3  Design  of  a-3D  wing:  The  third  example  considers  the  minimization  of  the  drag 
force  at  constant  lift  over  a  hypersonic  wing  (at  an  angle  of  attack  of  5°).  The  wing  is 
subject  to  a  minimum  thickness  constraint.  Figs.  3a  shows  the  Initial  geometry  and 
surface  mesh  on  the  wing,  Fig.  3b  shows  the  initial  and  final  wing  geometries  after 
80  design  cycles  (the  number  of  design  parameters  where  approximately  3575).  The 
design  was  mostly  dominated  by  the  minimum  thickness  contraint,  even  though  some 
small  twist  and  asymmetry  were  also  obtained.  The  objective  function  gradients  were 
computed  using  an  incomplete  continuous  adjoint  formulation.  The  design  process 
was  considered  finished  when  the  drag  and  the  lift  forces  did  not  have  any  further 
meaningful  variations.  The  drag  improvement  was  of  60%,  which  was  obtained  with  a 
lift  increment  of  9%. 


Fig.  3a:  Initial  geometry  and  surface  mesh  on  the  wing  (3575  nodal  points).  The  corner 
coordinate  values  are  also  shown 


Fig.  3b:  Initial  (dashed  line)  and  final  (continuous  line)  geometries.  A  small  twist  at 
the  final  geometry  can  be  noticed. 
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6.  CONCLUSIONS  AND  OUTLOOK 


We  developed  a  scheme  that  allows  for  the  direct  generation  of  clouds  of  arbitrary 
objects  with  the  same  degree  of  flexibility  as  advanced  unstructured  mesh  generators. 
We  implemented  adaptive  embedded  unstructured  grid  methods  in  our  flow  code.  This 
technique  has  been  extremely  successful  for  complex  fluid-structure  interaction  cases,  as 
well  as  some  external  aerodynamics  cases.  We  have  developed  a  new  numerical  method 
for  computing  three-dimensional  compressible  turbulent  flows  for  complex  geometries 
using  a  one-equation  turbulence  model  on  unstructured  grids.  The  method  uses  con¬ 
tainment  dual  control  volume  to  discretize  the  inviscid  fluxes  to  achieve  the  dramatic 
improvement  in  accuracy  on  a  highly  stretched  tetrahedral  mesh  and  the  matrix-free 
implicit  GMRES+LU-SGS  method  to  achieve  the  fast  convergency. 
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